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' Abstract 

in ' 

^vq . We develop a numerical solver for radiative transfer problems based on the weighted 

essentially nonoscillatory (WENO) scheme modified with anti-diffusive flux correc- 
tions, in order to solve the temperature and ionization profiles around a point source 
of photons in the reionization epoch. Algorithms for such simulation must be able to 
handle the following two features: 1. the sharp profiles of ionization and temperature 
at the ionizing front (I- front) and the heating front (T- front), and 2. the fraction of 
neutral hydrogen within the ionized sphere is extremely small due to the stiffness 

i 

of the rate equations of atom processes. The WENO scheme can properly handle 
these two features, as it has been shown to have high order of accuracy and good 
convergence in capturing discontinuities and complicated structures in fluid as well 
as to be significantly superior over piecewise smooth solutions containing disconti- 
nuities. With this algorithm, we show the time-dependence of the preheated shell 
around a UV photon source. In the first stage the I- front and T-front are coincident, 
and propagate with almost the speed of light. In later stage, when the frequency 
spectrum of UV photons is hardened, the speeds of propagation of the ionizing and 
heating fronts are both significantly less than the speed of light, and the heating 
front is always beyond the ionizing front. In the spherical shell between the I- and 
T-fronts, the IGM is heated, while atoms keep almost neutral. The time scale of 
the preheated shell evolution is dependent on the intensity of the photon source. 
We also find that the details of the pre-heated shell and the distribution of neutral 
hydrogen remained in the ionized sphere are actually sensitive to the parameters 
used. The WENO algorithm can provide stable and robust solutions to study these 
details. 
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shock waves 
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1 Introduction 



The profile of HII region around an isolated source of UV photons is an old 
topic in astrophysics. A classical result was given by Stromgren (1939), who 
showed that the profile of the spherical HII region of a point source embed- 
ded in a uniformly distributed hydrogen gas with number density n at radial 
coordinate r = can be approximately described by a step function as 

, 1, if r < R s 

Mr)~B[R a -r\ = { (1) 

0, if r > R s 

where the fraction of ionized hydrogen fmiij') = ^Hii( r )/^, ^Hii( r ) being the 
number density of ionized hydrogen. That is, hydrogen gas is sharply divided 
into two regions: within a sphere with Stromgren radius R s around the source, 
hydrogen is fully ionized, while outside the sphere hydrogen atoms remain 
neutral. R s is determined by the balance between ionization and recombination 
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where N is the emission of ionizing photons of the source, and cub is the recom- 
bination coefficient of HII (Osterbrock & Ferland 2005). The sharp boundary 
at r = R s is the ionization front (I-front) separating the HII and HI regions. 

The problem of the Stromgren sphere has once again attracted many studies 
recently, because the formation of HII regions around high redshift quasars, 
galaxies and first stars is crucial to understanding the evolution of the reion- 
ization (Cen & Haiman 2000; Madau & Rees 2000; Ciardi et al. 2001; Ricotti 
et al. 2002; Wyithe & Loeb, 2004; Kitayama, et al. 2004; Whalen et al. 2004; 
Yu, 2005; Yu & Lu 2005; Alvarez et al. 2006, Iliev et al. 2006). Unlike the 
static solution eq.(l), new studies focus on the dynamical behavior of the ion- 
ized sphere, such as the time-dependence of the ionization profile fmi{t,r), 
the propagation of the I-front. Besides the HII region and the I-front, a high 
kinetic temperature region also exists around the UV photon source due to 
the photon- heating of gas. Similar to the I-front, there is also a T-front sep- 
arating heated and un-heated gas. The temperature profile T(t,r) and the 
T-front are important for probing reionization. For instance, the region with 
high kinetic temperature T and low /fin would be the 21 cm emission region 
associated with sources at the reionization epoch (Tozzi et al. 2000; Wyithe 
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et al. 2005; Cen, 2006; Chen & Miralda-Escude 2006). Although the possi- 
ble existence of a 21 cm emission shell around high redshift quasars and first 
stars has been addressed qualitatively or semi-analytically in these references, 
a serious calculation seems to be still lacking. 

Many numerical solvers for the radiative transfer equation have been proposed 
(Ciardi et al. 2001, Gnedin & Abel 2001, Sokasian et al. 2001, Nakamoto et al. 
2001; Razoumov et al. 2002, Cen 2002, Maselli et al. 2003, Shapiro et al., 2004; 
Rijkhorst et al. 2005; Mellema et al. 2006; Susa 2006, Whalen & Norman 2006). 
These solvers provide numerical results of the I- front. However, the results are 
still diverse due to the usage of different approximations. Some of the results 
show that the time scale of the I-front evolution is sensitively dependent on the 
intensity of source (White et al. 2003), while some yield intensity-independent 
evolution (e.g. Mellema et al. 2006). This is because the retardation of photon 
propagation is ignored in the later, while the recombination is ignored in the 
former. Therefore, it is worth to re-calculate this problem without above- 
mentioned assumptions. It has been pointed out that to study the dynamical 
features of the I- and T-fronts, one would need to apply high-resolution shock- 
capturing schemes similar to those developed in fluid dynamics (Razoumov & 
Scott 1999). The finite difference WENO scheme is an algorithm satisfying 
this requirement. 

Moreover, although the fraction of the remained neutral hydrogen within the 
ionized sphere is extremely small, it is not zero. The small fraction is important 
to estimate the Ly-a photon leaking at high redshifts. Therefore, a proper 
algorithm for the dynamical properties of the ionized sphere should be able to, 
on the one hand, effectively capture sharp profile of ionization and temperature 
around the I- and T-fronts, and, on the other hand, give a precise value of 
the remained neutral hydrogen between the discontinuities. This can also be 
satisfied by the WENO algorithm. 

The WENO algorithm has proved to have high order of accuracy and good 
convergence in capturing discontinuities and complicated structures in fluid 
as well as to be significantly superior over piecewise smooth solutions con- 
taining discontinuities (Shu 2003). We have shown that the WENO algorithm 
is indeed effective for solving radiative transfer problem with discontinuities 
with high accuracy. For instance, it can follow the propagation of a sharp I- 
front and the step function cut-off of retardation (Qiu et al. 2006). We now 
develop this method to solve both ionization and temperature profiles. It is 
not a trivial generalization. Because the rate equations of heating-cooling and 
ionization-recombination are stiff, the time integration of the WENO scheme 
cannot be directly implemented. A proper strategy to save computational cost 
on time integration will be developed. A long-term motivation of this work, 
as we mentioned in (Qiu et al. 2006), is to develop a WENO solver of hydro- 
dynamic/radiative transfer problems, similar to the development of a hybrid 
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algorithm of hydrodynamic/N-body simulation based on WENO scheme (Feng 
et al. 2004). 

The paper is organized as follows. Section 2 describes the problem and equa- 
tions needed to be solved. Section 3 presents the WENO numerical scheme. 
Section 4 gives the solutions of the temperature and ionization profiles and 
the evolution of energy spectrum of photons. A discussion and conclusion are 
given in Section 5. Details of the atomic processes are listed in the Appendix. 



2 Basic Equations 



To demonstrate the algorithm, we consider the ionization of a uniformly dis- 
tributed hydrogen gas in space x with number density n by a point UV photon 
source located at the center |x| = r = 0. Adding helium component in the gas 
is straightforward and will not change the algorithm for radiative transfer. If 
the time scale t of the growth of the ionized sphere is less than 1/H(t), H(t) 
being the Hubble parameter, the expansion of the universe can be ignored. 
The radiative transfer equation of the specific intensity J{t, r, u) is then (see 
the Appendix) 

where v is the frequency of photon. The source term, S, is given by 

S(t,\^M = E(u)5(x) (4) 

where E{y)dv is the energy distribution of photons emitted by the central 
source per unit time within the frequency range from v to v + dv. We assume 
the energy spectrum of UV photons to be of a power law E[y) = Eo(h>o/v) a , 
and uq is the ionization energy of the ground state of hydrogen hu Q = 13.6 
eV. Integration of E over v gives the total intensity (energy per unit time) of 
ionizing photons emitted by the source, E = E{v)dv = E u /(a — 1). 

The absorption coefficient of eq.(3) is k u = <r(z/)nHi(t, x), where the cross 
section a(u) ~ a (is /is) 3 and a = 6.3 x 10 -18 cm 2 . The evolution of the 
number density of neutral hydrogen HI, nm(t, r), is governed by the ionization 
equation, 

~~rr~ = a an n efmi — r 7 m/Hi — r c Hin e /Hi (5) 

where /hi = n}n(t,r)/n is the fraction of neutral hydrogen, and n e is the 
number density of electrons. Obviously, fm(t,r) + /Hii(^ r ) = 1- In eq.(5), 
ami is the recombination coefficient and r e ni is the collision ionization rate. 
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The photoionization rate T 7 Hi(t, r) is given by 



The kinetic temperature of baryon gas is determined by the equation 

jrri 

nk B -- = H-n 2 C (7) 
dt w 

where ks = 1.38 x 10~ 16 erg K~ l is the Boltzmann constant and the temper- 
ature T is in unit of K. The details of the heating H and cooling C are given 
in the Appendix. 



3 Numerical Algorithm 



3.1 Dimensionless variables 



In the numerical implementation, it is convenient to introduce the dimension- 
less variables of time, space and frequency defined by t' = caoiit, r' = a^nr 
and v' = u/uq. I/ctqU is the optical depth of ionizing photon in neutral hy- 
drogen gas with density n. Therefore, t' and r' are respectively, the time and 
distance in units of mean free flight time and mean free path of ionizing pho- 
ton his® in the non-ionized background hydrogen gas n. For the ACDM model, 
n = 1.88 x 1(T 7 (1 + zf cm" 3 , where z is redshift, t' = 0.89(1 + z)- 3 t Myrs 
and r' = 0.27(1 + z)~ 3 r Mpc. Correspondingly, the intensity is rescaled by 
Jdv = {l/r 2 ){hv /aln)J'dv'. Thus, eqs.(3), (5), (6) and (7) become 

and 

Ca °~dY = aHII ^n ~ ~^f^ HI ~~ r oHi( 1 ~~ /hi) /hi- ( 9 ) 
r 7ffl (t,r) _ 1 r ,J'(t, r ,i/) 1 



n 



dt 1 n 

where we have assumed n e /n = fun and 

1 hu Q £ f°° v' - 1 
H 



r' 2 J i v v' 6 

C ^B§ = ^H-C. (11) 



n 2 r' 2 



hi l — J>d »' (12) 
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The point source condition eq.(4) requires E{y)dv = 47r(/iz/ /n<7o) J'(t, 0, v'^dv' . 
If we take J'(t, 0, z/) = J (l/V') a , the total intensity of the source is given by 

E = , 4nhU ° 3 Jo = 5.8 x 10 42 ^- (4-) erg/s. (13) 

{a-l)nal a-l\l + z r J VIO- 3 / &/ v ; 

In our numerical calculation, we solve the system of equations (8), (9) and 
(11) for the specific intensity J', the fraction of the neutral hydrogen / H i and 
the temperature T as functions of the radius r', frequency v' and time t' . We 
will drop the prime in the variables J', r', v' and t' in this section below, 
when there is no ambiguity, and keep prime in the variables in showing the 
numerical results. 

To solve the radiative transfer equation, we adopt the fifth-order finite differ- 
ence WENO scheme with anti-diffusive flux corrections. The fifth-order finite 
difference WENO scheme was designed in (Jiang & Shu 1996) and the anti- 
diffusive flux corrections to the high order finite difference scheme was designed 
in (Xu & Shu 2005). The objective of the anti-diffusive flux corrections is to 
sharpen the contact discontinuities in the numerical solution of the WENO 
scheme as well as to maintain high order accuracy. A fourth order quadrature 
formula is used in the computation of integration in equations (10) and (12). 
The third order TVD Runge-Kutta time discretization is used in time inte- 
gration for the system of equations (8), (9) and (11). We now describe our 
numerical algorithm in more detail. 

3.2 The computational domain 

The computational domain is (r, v) e [0,r maa; ] x [l,z/ ma J, where r max and 
Vmax are chosen such that J(r, z/, t) ~ for r > r max or v > u, max . In our 
computation, r max = 1200 and u max = 10 6 . The computational domain is 
discretized into a uniform mesh in the r-direction and into a smooth non- 
uniform mesh in the //-direction. The uniform mesh in the r-direction is 

r { = iAr with Ar = r max /N r , i = 0, N r 

and the non-uniform mesh in the //-direction is taken to be 

Vj = & with = j'A£, A£ = \og 2 u max /N u , j = 0, N v 

which is allowed because only integration, e.g. in equations (10) and (12), is 
involved in the computation with respect to the v- variable. In our model, large 
v contributes little in equations (10) and (12), therefore the non-uniform mesh 
is designed in a way such that the mesh becomes coarser for larger v. 
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3.3 Approximation to the spatial derivative 



To approximate the spatial derivative in equation (8), the WENO scheme 
with ant i- diffusive flux corrections is used. Specifically, to calculate dJ/dr, 
the variable v is fixed and the approximation is performed along the r-line 

|: J (^^^)^^(% + i/2-^-i/ 2 ) (14) 

where the numerical flux h? +1 , 2 is obtained with the procedure given below. 
We can use the upwind fluxes without flux splitting in the fifth order WENO 
approximation because the wind direction is fixed (positive). To obtain the 
sharp resolution of the contact discontinuities, the anti-diffusive flux correc- 
tions are used in our code. 

First, we denote 

hi = J{t n ,r u vj), i = -2,-l,...,N r + 2 

where n and j are fixed. The numerical flux from the regular WENO procedure 
is obtained by 

K+l/2 = W A ( +l/2 + ^2^+1/2 + 

where h^y 2 are the three third order fluxes on three different stencils given 
by 



f(i) 1, 7 11 

ll i+l/2 = g _ g + ~g 

£( 3 ) _ 1 h , 5 h l h 
n i+l/2 ~ 3 n i + Q i+1 ~~ Q i+2) 

and the nonlinear weights u m are given by 

& m _ 7/ 

z^=i (£ + pi) 

with the linear weights 7; given by 

1 _ 3 3 

7i — Jo' 72 ~ 5' 73 ~ To' 

and the smoothness indicators fli given by 
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13 1 

A = — (hi..* - 2h t ^ + li,) 2 + - (hi..* - 4^_x + 3k,) 2 

13 2 1 2 

^ 2 = 12 ~~ 2ki + hi+1 ^ + 4 _ 

13 2 1 2 

^ 3 = 12 ^ ~~ 2ki+1 + hi+2> + 4 ^ ~~ 4/l * +1 + ^ +2 ^ ' 

e is a parameter to avoid the denominator to become and is taken as e = 1CT 5 
times the maximum magnitude of the initial condition J in the computation 
of this paper. The reconstruction of the finite difference WENO flux on the 
downwind side hf +1 , 2 is obtained in a mirror symmetric fashion with respect 

to Xi+1/2 as that for h~ +1 / 2 . 

The anti-diffusive flux corrections are based on the fluxes obtained from the 
regular WENO procedure. It is given by 



K+i/2 — h-i+i/2 + 

^minmod 



(15) 



hi — hi-i 



+ ^i-l/2 h i+1 / 2 , ll i+l/2 



'ht +1 /2 h: 



i+1/2 



where r\ = At/Ar is the CFL number and the minmod function is defined as 



minmod(a, b) 



0, if ab<0 

a, if ab > 0, | a |<| b 

b, if ab > 0, | b \<\ a \ 



(16) 



4>i in eq.(15) is the discontinuity indicator between and 1, defined as 



(pi 



Pi 



Pi + li 



where 



Oti 



+ 



li 



hi-i - hi | +C) 2 , 



with ( being a small positive number taken as 1CT 6 in our computation. u max 
and « min are the maximum and minimum values of hi for all grid points. With 
the definition above, we will have < (pi < 1. (pi — 0(Ar 2 ) in the smooth 
regions and (pi is close to 1 near strong discontinuities. The purpose of the anti- 
diffusive flux corrections is to improve the resolution of contact discontinuities 
without sacrificing accuracy and stability of the original WENO scheme. 
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3.4 High order numerical integration in v 



The integration in equations (10) and (12) is approximated by a fourth order 
quadrature formula 

/ f(x)dx = AxY,w J f(jAx) + 0(Ax 4 ) (17) 

3=30 



where u = jo Ax, and the weights Wj are given by 

3 7 23 

w io = g, w io+i = g> w h+2 = w jo+i = 1, for j > 2. 

Again we refer to (Qiu et al. 2006) for details of implementation. 



3.5 Time integration 



When considering the time integration for the system of equations (8), (9) and 
(11), we start with the third order TVD Runge-Kutta time discretization in 
(Shu & Osher 1988). For the system of ODEs u t = L(u) this time discretization 
reads 



u^ = u n + AtL(u n ,t n ) (18) 
M (2) = ^ M « + I( M (i) + AtL(uW)) (19) 

u n+l = -u n + -{u^ + AtL(u< 2 >)) (20) 



The Runge-Kutta method needs to be modified considering the modification 
on the anti-diffusive flux f a by 



J (1) = J n + AtL(J n ,t n ) (21) 
= r + ^AtL\J n ) + ^AtL(J^) (22) 

jn+l = J n + I AtL "(jn) + -AtL(J^) + -AtL(J^) (23) 

6 6 3 

where the spatial derivative dJ/dr term in operator L is defined by (14) with 
the anti-diffusive flux h given by eq.(15), and dJ/dr in the operator V is 
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defined by the modified anti-diffusive flux h as 

'4(/i< - 



"i+l/2 



\+i/2 + mmm °d 



r/ 



i— 1/2 



"i+l/2' "i+l/2 



i+l/2 



(24) 



h a 

"-•i+l/2) 



«/ bc> 0, | 6 |<| c 
otherwise 



and dJ/dr in the operator L" is defined by the modified anti-diffusive flux h a , 

'6(hi - hi-i) 



"i+l/2 



h 



i+l/2 



minmod 



r/ 



' ^i+1/2 ^i+l/2) ' 

i/ 6c> 0, 1 6 |<| c 
otherwise 



+^i-l/2 ^i+l/2' ""i+l/2 



"i+l/2> 



(25) 



Here b = (hi - h^/r] + ^_ 1/2 - fe i+1/2 , c = fe+ 1/2 - ^ +1/2 . 

The difficulty of a direct implementation of the scheme lies in the stiffness of 
equations (9) and (11). Especially for strong sources, when J(r=0) is large, 
one needs very small time step At to guarantee the stability of the numerical 
scheme, therefore huge computational cost for long time integration. 

By observing that 

• Though (9) and (11) are stiff ODEs, (8) involving the spatial derivative is 
not a stiff equation. 

• The WENO procedure with the ant i- diffusive flux corrections in approxi- 
mating the spatial derivative is the major cost at each time step. 

• Implicit numerical method in time evolution has milder time step restriction 
than the explicit one. 

we have settled down with the following strategies to save the computational 

cost 



• We use different time scales to solve J in (8) and to solve fni and T in 
(9) and (11). The time step for J, say At, is larger, while the time step for 
fni and T, say St = ^ with N t being the number of small time steps in a 
large time step, is much smaller. By doing this, we take advantage of the 
non-stiffness of (8) and eliminate the cost of the WENO procedure with flux 
corrections for each small time step. 

• In each small time step St, we use a semi-implicit numerical method for 
strong sources to proceed in time for fni and T, which greatly releases the 
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severe time step restriction for the sake of the stability of the scheme. The 
number of small time steps St in one large time step At is thus greatly 
reduced, hence the saving of the computational cost. 

More precisely, we modify the third order TVD Runge-Kutta time discretiza- 
tion in the following procedure 

(1) We have the initial condition of J, f HI and T at t=0. Let us denote 

fm = (fmi)N r xi = (f HI (xi,t = 0)) NrXl ; T° = (Tf) NrXl = (T(x t ,t = 0)) NrXl . 

(2) For a large time step At, we evolve J n , f HI and T n by three inner stages 
to J n+1 , and T n+1 . 

• J n+1 : 

J (i) j j(2) ; jn+i are U pd a ted according to (21), (22), (23). 

fn+l rpn+1. 

• J HI 1 1 

Let us denote 

h = (h 1 ,h 2 ) = (f HI ,T) 

and 

L(J, f HI , T) = (Lx( J, f HI , T), L 2 ( J, f HI , T)) 

= C J_{ rhs of (9)}, — { r/w of (11)}) 
\C(7o c<JqKb / 

Instead of (18), f m and T^> are evolved by 

where h(N t ) is the evolution of the solution by the Euler-forward method 
for N t times, with the Euler-forward method at each time step defined 
as 

h {]+1) = h i3) + 5tl{J\h {]) )- (26) 
Instead of (19), the second inner stage h^—(f m ,T^) is updated by 

- h m = h n + \m 

4 4 

with h^ 2 '^ = h(N t ) being the evolution of the solution by the Euler- 
forward method for N t times based on as 

kj+i) = hj) + 6tL (J (1 \h U )^ ( 27 ) 
Finally, instead of (20), f H + j l and T n+1 are computed as 

- h n+l = hn + 2^(3') 
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with M 3 ') = /i(jv t ) obtained in a similar manner as in the previous step 
but based on 



h U+V = h U) + 5tL ( j(2) ' h ti))- 



(28) 



(3) When the source is strong, (9) and (11) suffer from severe time step 
restriction, we use a semi-implicit scheme instead of (26) in our code to 
release the severe time step restriction and to save computational cost. 
The procedure of the implementation is the same as described above 
except that (26) is replaced by 



(29) is computed by solving a quadratic equation with the root located 
between and 1, since L\ is a quadratic function of fni, and (30) is 
computed by the Newton iteration method. (27) and (28) are modified 
in a similar way as (26). 

3.6 Boundary conditions 

The boundary conditions are implemented as follows 

• Inflow boundary condition at r=0: 

Ji,j = Jo,j, for i = 0,-1,-2. 

• The boundary condition at r = r max : 

JN r +i,j = Jn t -ij, f or i = 0, 1, 2 

3. 7 Convergence study of the numerical scheme 

In this subsection, we perform a grid refinement convergence study for the 
numerical scheme we proposed above, to assess the accuracy of the computa- 
tional result for the typical mesh sizes that we will use in next section. 

We test our numerical schemes for the case with the strongest source intensity, 
i.e. E = 5.8 x 10 45 erg s" 1 . This is the toughest situation for our numerical 
simulation, as the scheme suffers severe time step restriction due to the stiffness 
of the ODE. We use both the multi-time-scale strategy and semi-implicit time 
discretization described in the previous subsections, and compare the results 
obtained with 4000 points in r and 200 points in is, which is the typical mesh 
size used in next section, and with a coarser mesh consisting of 2400 points 
in r and 100 points in u, in Figure 1. The results with these two mesh sizes 



fHi(j+i) = fm(j) + StL^J" 1 , fHi(j+i),T(j)) 
T (i+i) = T (j) + StL 2 (J n , f HI Q +1 ),T {j+1) ) 



(29) 
(30) 
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Fig. 1. The profiles of f m (t, r') (left) and T(t, r') (right) with the source of 5.8 x 10 45 
erg s _1 at time t = 0.89 Myrs (t' = 1000). The power-law frequency spectrum has 
the index a = 2, and the reionization redshift is taken to be 1 + z = 10. The solid 
line is the numerical result with the typical mesh N r = 4000, N„ = 200, and the 
circles indicate the numerical result with a coarser mesh N r = 2400, N v = 100. 




Fig. 2. The profiles of f m (t, r') (left) and T(t, r') (right) with the source of 5.8 x 10 39 
erg s _1 at time t = 0.09 Myrs (t 1 = 100). The power-law frequency spectrum has 
the index a = 2, and the redshift is taken to be 1 + z = 10. The numerical mesh 
is N r = 2400 and N v = 200. The solid line is the numerical solution with single, 
extremely small time step, and the circles indicate the numerical solution with the 
multi-time-scale strategy. 

match very well, indicating that our numerical results are already numerically 
convergent at the coarser mesh. The numerical results reported in next section 
with the finer mesh should therefore be reliable. 

We also test the accuracy of our multi-time-scale strategy by comparing the 
numerical results of evolving (8), (9) and (11) together with single, extremely 
small time step, and that of updating (8), (9) and (11) by the multi-time-scale 
strategy with (5t, At). These two approaches produce numerical solutions in 
good agreement, see Figure 2. 

Additional grid refinement study, not reported here to save space, has been 
performed to assure the numerical convergence of the results reported in next 
section. 
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4 Results 



The results in this section are obtained with the numerical algorithm described 
in the previous section. The code is stable for all the cases reported in this 
section. As indicated in the previous section 3.7, and we have performed grid 
refinement study for some representative examples to ensure that the results 
reported are numerically converged solutions. 

In our calculation, we use the explicit time discretization to solve (9) and (11) 
when E = 5.8 x 10 39 erg s" 1 , and semi-implicit scheme when E = 5.8 x 10 41 
erg s -1 , 5.8 x 10 43 erg s _1 and 5.8 x 10 45 erg s _1 in each of the small time step 
St in the multi-time-scale strategy. 

4-1 Profiles of fui(t,r) andT(t,r) 

We calculate the profiles of T(t,r) and fm(t, r ) around a point source. The 
result is presented in Figure 3, which shows the time-dependence of the profiles 
fm{t, r) and T(t, r) for sources with intensity, respectively from top to bottom 
panels, E = 5.8 x 10 39 , 5.8 x 10 41 , 5.8 x 10 43 and5.8xl0 45 erg s" 1 with a = 2 and 
1 + z — 10. The emission rate of the number of photons is roughly N ~ 10 50 , 



Figure 3 shows that the ionized sphere generally has a sharp I- front. The time- 
dependent profile fm(t, r) may still be approximately described as a Stromgren 
sphere with time-dependent size R(t) as 



The code can well reveal the propagation of the I-front r = R(t). Unlike 
the solution eq.(l), fm(t,r) is not zero in r < R(t). Although the neutral 
hydrogen fm(t,r) remained in r < R(t) is very small, it cannot be ignored, 
because the heating of gas is given by the hard photons absorption of neutral 
hydrogen HI [eq.(A8)]. The high temperature within the ionized sphere is 
actually maintained by the ionization of the small fraction therein by all the 
ionizing photons [eq.(31)]. Figure 3 indicates that fm(t,r) within r < R(t) 
is sensitive to the intensity of the source. This is important for investigating 
leakage of the Ly-a photon. 

Figure 3 shows that the temperature profiles T(t, r) within the I-front r < 
R(t) are simple. The temperature keeps around a constant 10 4 — 10 5 K. An 
interesting feature is that for strong sources, the temperature at the center 



10 52 , 10 54 and 10 56 s" 1 . 




(31) 
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Fig. 3. The profiles of fm(t,r) (left) and T(t,r) (right), in which r and r' are the 
physical and dimenssionless distance, respectively. From top to bottom: for sources 
of E = 5.8 x 10 39 , 5.8 x 10 41 , 5.8 x 10 43 and 5.8 x 10 45 erg s" 1 at time t = 0.18 
(dash dot dot), 0.36 (long dash), 0.53 (dash dot), 0.71 (dash), and 0.89 (solid line) 
Myrs. The power-law frequency spectrum has the index a = 2, and the redshift is 
taken to be 1 + z = 10. The mesh has TV r =4000, iV„=200 grid points. 
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(r ~ 0) is not the highest and is even lower than the constant temperature. 
This is because the ionizing (soft) photon flux at the center is very strong, the 
number of /m(r ~ 0) is extremely low (see Figure 3), and then, the heating 
rate is lower. 

A common feature of T(t, r) is that the size of heated region is generally 
larger than that of ionized region at all time. That is, hydrogen gas is already 
significantly heated before being ionized. The gas is firstly heated up, and then 
ionized. The region between the T-front and the I-front is a pre-heating layer, 
in which the temperature can be as high as T ~ 10 3 K, or even larger. For a 
source of E = 10 43 erg s _1 at (1 + z) = 10, the physical size of the ionized 
region at time 0.36 Mys is about 0.06 Mpc, while the physical size of the 
region with temperature larger than 10 3 K is about 0.08 Mpc. When t = 0.89 
Mys, the two physical sizes are, respectively, about 0.1 Mpc, and 0.16 Mpc, 
or comoving size 1 Mpc, and 1.6 Mpc. Therefore, for strong sources, the time 
scales and comoving length scales of the pre-heating layer are in the range of 
cosmo logical interest. 

An important feature of the ionization profile is that fm(t,r) can be ap- 
proximately written as fm(t,r) = fm(r)9[r — R(t)], where fm(r) is time- 
independent, and R(t) is only a function of t. That is, all the time-dependence 
of fm(t,r) can be described by the I-front R(t). In the region r < R(t), the 
ionization profile fm(t,r) = fm( r ) is time-independent. It can be given by a 
static or stationary solution of the radiative transfer equation. This property 
has been applied in some codes for radiative transfer, in which the ionization 
field is given by the static solution of the radiative transfer equation for a 
given matter distribution. 

For the temperature profile, one might also approximately define a T-front 
function by a step function as 8[t — R(t,T)]. However, unlike the I-front, the 
temperature profile can not be rewritten as T(t,r) = T(r)9[t — R(t,T)], where 
T(r) is time-independent. That is, the r-dependence of the function T(t, r) can 
not be separated with t. Especially, the function T(t,r) in the range between 
the I- and T-fronts actually always depends strongly on t. 



4-2 Speed of the I- and T-fronts 

Figure 3 shows that, for a very strong source E = 10 45 erg s _1 , we have 
R(t) ~ ct, i.e. the ionizing front moves with a speed close to the speed of 
light. For a weak source E = 10 38 erg s _1 , the profile fm{t, r) does not show a 
significant expansion of the I-front if time t is larger than a few of free flight 
times of ionizing photons. This result is also evident in Figure 4 on r 90 (t) 
versus t. r 90 (t) is defined by /Hi(^ r 9o) — 0.90. It denotes the size, within 
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Fig. 4. left panel: the evolution of rgo(t), which is the solution of /Hi(^ r 9o) = 0.90. 
The light solid line is rgo = ct. Right panel: drgo/dt vs. t. Other parameters and the 
mesh in the numerical simulation are the same as those in Figure 3. The sources are 
taken to be E = 5.8 x 10 39 (dash dot dot), 5.8 x 10 41 (dash dot), 5.8 x 10 43 (dash), 
and 5.8 x 10 45 (solid line) erg s^ 1 . Time t' is dimensionless. 

which, i.e. r < rgo, 90% hydrogen are ionized. Therefore, it can be used for the 
I-front. For a strong source, such as E = 10 45 ergs s _1 , we have approximately 
r 90 = ct, which implies that the ionizing region grows with an ionizing front 
propagating with almost the speed of light. For weak sources, r 90 (t) are also 
following t at very small t, but become r 90 (t) < ct when t is large. The weaker 
the sources, the earlier the stage of rgo(t) < ct takes place. This point is more 
clearly shown in the right panel of Figure 4. We see that the speed dr$o(t)/dt is 
close to c when t is small, and then the speed decreases with t by a power law 
dr 90 (t)/dt oc t~ p with (3 > 2/3. One can define a time U, larger than which, 
the speed dr$o(t)/dt starts to decrease with t by the power law. At time t > t it 
the ionizing sphere is still expanding, but very slowly. It will finally approach 
a solution, of which the ionization equilibrium is approximately established, 
and the ionized sphere becomes time- independent. The time ti depends on the 
source intensity. The stronger the source, the larger the time tj. 

Obviously, the speed of the propagation of T-front cannot be larger than the 
speed of light, and therefore, the T-front will approximately coincide with the 
I-front when t < t^. Only in this period, the ionized sphere is the same as 
the heated sphere. When t > U, the T-front starts to exceed the I-front, and 
the pre-heating layer is formed. Therefore, the formation of pre-heated layer 
happens later for a stronger source. 

Figure 5 gives a comparison of fm(t,r) and T(t,r) at time t = 0.89 Myrs for 
sources with different intensity. We can see that the pre-heating layer has been 
well established at time t = 0.89 Myrs for all sources with E < 10 43 erg s -1 , 
while the T-front of the source with E = 10 45 erg s _1 is still about the same as 
the I-front. One can expect that for a source of E > 10 45 erg s -1 , a preheated 
layer will be formed at time t > 0.89 Myrs and on comoving distance r > 2.7 
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Fig. 5. Profiles of /hi(^) (left panel) and T(r) (right panel) at time t = 0.89 Myrs 
for sources with intensity E = 5.8 x 10 39 (dash dot dot), 5.8 x 10 41 (dash dot), 
5.8 x 10 43 (dash), and 5.8 x 10 45 (solid line) erg s _1 . r and r' are the physical and 
dimensionless distance, respectively. The power-law frequency spectrum has index 
a = 2. The redshift is taken to be 1 + z = 10. The mesh in the numerical simulation 
is the same as that in Figure 3. 

Mpc from the source. 



4-3 Spectral hardening 

The formation of the pre-heated layer with high T and high / HI is due to the 
lack of soft photons beyond the I-front, while hard photons are still abun- 
dant in that region. In other words, the energy spectrum of the photons is 
significantly hardened around the I-front. We now directly demonstrate the 
evolution of the photon frequency spectrum, as our code can effectively reveal 
the evolution of radiations in the z/-space. 

The left panel of Figure 6 gives the frequency spectra of 1.) source E = 
5.8 x 10 43 erg s^ 1 at time t = 0.89 Myrs and physical distance r = 0.08, 
0.11, and 0.14 Mpc, and 2.) source E = 5.8 x 10 45 erg s" 1 at time t = 0.89 
Myrs and physical distance r = 0.19 and 0.24 Mpc. We can see a significant 
r-dependence of the frequency spectrum. At a small r the spectra are still 
almost the same as the original power-law spectrum u~ a with a — 2, while at 
large r, they significantly deviate from the original power-law. All photons of 
v < 10uq, are exhausted within r < 0.24 Mpc. At high frequency v > 50u , 
the spectra are still of power-law with a — 2, while at v < 50u they are 
substantially dropped. It shows a peak at 1 < u/u < 10, and looks like a 
spectrum with self-absorption. However, unlike a self-absorption spectrum, 
the position of the peak is not fixed in the frequency-space, it moves to higher 
frequency with time. Namely, the photon energy spectrum is harder at later 
time. 



18 




10 4 (eV) 



10 4 (eV) 




Fig. 6. J'(t,r,u) vs. logv/vo at time t = 0.89 Myrs. 1.) Left panel: for the source 
E = 5.8 x 10 43 erg s _1 , at physical distance r = 0.08 (dash dot), 0.11 (dash) and 
0.14 (solid line) Mpc; and 2.) Right panel: for the source E = 5.8 x 10 45 erg s _1 at 
physical distance r = 0.19 (dash) and 0.24 (solid line) Mpc. Redshift is 1 + z = 10. 
The mesh in the numerical simulation is the same as that in Figure 3. 
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Fig. 7. a vs. logz//Vo for the source E = 5.8 x 10 45 erg s _1 at physical distance 
r = 0.19 (dash) and 0.24 (solid line) Mpc. Redshift is 1 + z = 10. The mesh in the 
numerical simulation is the same as that in Figure 3. 



The spectral hardening can be measured by the index of power law denned as 

9 In J 



a = 



dlnu 



(32) 



Figure 7 plots a vs. v for the frequency spectra of Figure 6. Obviously, in the 
band vjv§ < 10, a becomes smaller at larger r. Both Figures 6 and 7 show 
that the frequency spectra of photons depend strongly on t and r when t and 
r are close to or inside the pre-heating layer. 
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5 Discussions and conclusions 



We have described WENO scheme which is able to solve the phase-space 
distribution function of photons in an isolated ionized patch around individual 
source in the early stage of reionization. This algorithm can produce robust 
results for the propagation of the I- and T-fronts. It can also give stable results 
for the time- dependent distribution of the small fraction of neutral hydrogen 
within the I-front. We have developed the method to deal with the stiffness of 
rate equations in time integration. Consequently, the computational speed is 
acceptable. This algorithm can be applied to problems with a wide range of 
intensity of sources, from 10 39 to 10 45 erg s -1 with a power law spectrum at 
redshift 1 + z — 10. Since this algorithm treat the frequency space as well as 
physical space, it can also be used for sources with other frequency spectrum. 
We have not considered helium and secondary effects in the atomic processes 
yet, however the algorithm has no practical difficulty to include these factors. 

We show that a common feature of the UV photon sources in the reionization 
epoch is to form a preheating region. In the first stage the I- and T-fronts in 
baryon gas are coincident, and propagate with the speed about the same as the 
speed of light. When the frequency spectrum of the UV photons is hardened, 
the evolution enters into the second stage. The propagation speeds of both 
the ionizing and heating fronts are less than that of light, but the T-front is 
always moving faster than the ionizing front. In the spherical shell between the 
I- and T-fronts, the kinetic temperature of gas can be as large as T ~ 10 3 " 4 , 
while atoms are almost neural. Obviously, the shell would be of interest in 
the search for the 21 cm emission from the reionization epoch. The details 
of the preheating region are sensitively dependent on the parameters of the 
problems. For instance, the radius and thickness of the preheated shell do not 
show a scaling relation on their dependence of the source intensity. Therefore, 
an algorithm, which can properly handle models with various parameters is 
necessary. 

Using these results, one can make comments on some numerical solvers used 
for the radiative transfer equations. Several solvers are based on the approxi- 
mation of omitting the time-derivative term of the radiative transfer equation 
(Nakamoto et al. 2001; Cen 2002; Razoumov et al. 2002; Rijkhorst et al. 2005; 
Susa 2006). They calculate the static ionization field for a given uniform or 
non-uniform density fields of cosmic baryon matter. These codes essentially are 
generalizing the static solution eq.(l) to the case of inhomogeneous background 
distribution of baryon matter, and multiple sources. Such an approximation 
is useful to calculate the ionization field for each given inhomogeneous mass 
field of baryon matter, but it will no longer be a good approximation when the 
retardation effect due to photon propagation is important. For the problem 
of ionization profiles around a point source, the retardation of photon is not 
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always negligible. For instance, one finds the following parameters have been 
used in the model: the ionization and heating at comoving distance r = 10 h _1 
Mpc from a point source of age t = 0.5 x 10 6 yrs at redshift z = 30 (e.g. Cen, 
2006). These parameters already violated the retardation constraint r < ct. 

The retardation effect can be seen from the time-dependence of the I- front, 
r\(t). If the retardation effect is negligible, an analytical solution of ri(t) is 
given by e.g. Whalen & Norman (2006) 

r I (t) = R a [l-e- t ' t "°] 1 / 3 , (33) 

where R s = (37V / '4^0^11^) ^ 3 i s the static radius of the Stromgren sphere 
given by eq.(2), and t rec = I/ckhii^h * s recombination time scale. According 
to eq.(33), the time scale of the rj(t) evolution is t rec , which is independent 
of the source intensity E. That is, r\(t) ~ R s only if t ^> t rec , regardless 
N. However, Figure 4 shows that time scale of the r 90 (t) evolution is N- 
dependent. This result also holds if we replace r 90 (t) by, e.g. r 95 (t), and is true 
for a monochromatic source too. 

Considering the retardation effect, an analytical solution of ri(t) is approxi- 
mately given by the following algebraic equation (White et al. 2003) 



n(t) 



' (/•.(/)'•)' 

Ann/3N 



(34) 



According to eq.(34), when t <C t c — (3N /Aunc 3 ) 1 ^ 2 , the speed of the I-front, 
drj(t)/dt ~ c, and when t 3> t C) dri(t)/dt oc t -2 / 3 . It is qualitatively consistent 
with dr 90 (t)/dt shown in Figure 4. First, the time scale of the r 90 (t) evolution 
is approximately oc TV" 1 / 2 , the larger the N, the longer the t c . Second, dr\(t)/dt 
decreases with t by a power law drgo(t)/dt oc t^ 13 when t is large. However, the 
power law index (5 > 2/3. This is expected, because eq.(34) does not include 
the effect of recombination, which leads to a slower decrease than dr 90 (t)/dt. 

A common assumption used in eqs.(33) and (34) is that the distribution 
fm(t,r) is described by a step function like eq.(l). Obviously, it ignores the 
neutral hydrogen HI probably remained within ri(t). The tiny fraction of HI 
may also be hardly calculated well by Monte-Carlo codes (Ciardi et al. 2001; 
Maselli et al. 2003), which yield large numerical errors due to Possion shot 
noise. 



The WENO algorithm has been successfully applied to kinetic equations of 
the distribution function in the phase space with one or two spatial dimensions 
and two or three phase space dimensions with acceptable computational speed 
(Carrillo et al. 2006). We believe that it is not difficult to implement the 
WENO algorithm for radiative transfer problems with similar dimensions in 
the phase space. 



21 



In our calculation, the evolution the cosmic baryon gas is not tracked by 
the hydrodynamic equations, but is simply assumed to have a uniform dis- 
tribution with density n. This treatment would be reasonable if the typical 
time scale of the relevant hydrodynamic effects are less than that of the I- 
and T-fronts. For instance, dynamical effects associated with sonic propaga- 
tion would be negligible in solving the propagation of the I- and T-fronts. 
Of course, one can expect that richer results will be yielded if the WENO 
scheme for radiative transfer problems can be incorporated with the Euler hy- 
drodynamics. Since the WENO scheme for the cosmological hydrodynamical 
simulation has already been well established, it would be possible to develop 
a unified radiation/N-body/hydrodynamics code for cosmological problems. 
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A Equations of Radiative Transfer 

For an ionized sphere associated with a point photon source, the radiation 
transfer (RT) equation is (Bernstein, 1988, Qiu et al. 2006) 

% +1 al-^- H (^- 3J )- k " J + S - < A1 » 

where J(i, x, z/, n,j) is the specific intensity, a the cosmic factor, H = a/ a, v 
the frequency of photon and n,i a unit vector in the direction of photon prop- 
agation. In eq.(A.l), we take c = 1. k u and S are, respectively, the absorption 
and sources of photons. The absorption coefficient of eq.(A.l) is given by 

k v = <r(v)n m (t, x) (A.2) 

where the cross section a(v) = 6.3 x 10~ 18 (v /is) 3 cm 2 . 

The number density of neutral hydrogen HI, nm(t, x), is determined by 

= ttHII^e/ffll — r 7 Hl/ffl — r e Hl7l e /HI, (A. 3) 

Relevant parameters are taken from Theuns et al. (1998) as 
1. recombination coefficient 

a m i = 6.30 x 1Q- U T- 1/2 T 3 - 2 /(l + T 6 a7 ), (A.4) 
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where T is temperature, and T n = T/10 n . 

2. collision ionization 

r eHI = 1.17 x i -iOT 1 / 2 e- 157809 - 1 / T (l + Tl 12 )- 1 (A.5) 

3. photoionization 

r 7 m(t, r) = r dv J{t ^ V) a{v). (A.6) 

The temperature T is determined by the equation 

nk B — = H- n 2 C. (A.7) 
at 

4. heating rate 

H = n HI dvJ{t,r,v)a{v) (A.8) 

where Ylvq = 2.176 x 10 -11 ergs. 

5. cooling. Since only the recombination cooling is important, we have 

C = 8.70 x 10- 27 T 1 / 2 r 3 - - 2 (l + T 6 a7 )- 1 [l - / ffl ] 2 (A.9) 
+ 1.42 x 10" 27 T 1/2 [1 - /hi] 2 

+ 2.45 x i -2i T i/2 e -i57809.i/r (1 +t 5 1/2 )- 1 (1 - /hi)/hi 
+ 7.5 x 10- 19 e- 118348 / r (l + T 5 1/2 )- 1 (l - /hi)/hi 

where T n = T/10™. The terms on the r.h.s. of eq.(A.9) are, respectively, the re- 
combination cooling, collisional ionization cooling, collisional excitation cool- 
ing and bremsstrahlung. Both H and C are in the unit of ergs cm 3 s _1 . 
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